Motion of dark solitons in trapped Bose-Einstein condensates 
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We use a multiple time scale boundary layer theory to derive the equation of motion for a 
dark (or 'grey') soliton propagating through an effectively one-dimensional cloud of Bose-Einstein 
condensate, assuming only that the background density and velocity vary slowly on the soliton 
scale. We show that solitons can exhibit viscous or radiative acceleration (anti-damping), which we 
estimate as slow but observable on experimental time scales. 
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The success of the Gross-Pitaevski mean field theory in describing experimentally observed dilute Bose condensates 
shows that one really can persuade a large number of particles to behave as a field. There is thus a pleasant 
, circularity in investigating situations where this field in turn behaves in a particle-like manner, in that it contains 
PsJ ' topological defects or solitons. fn this paper we discuss one particular particle-like configuration of the Gross-Pitaevski 
^ ' mean field, namely the one-dimensional dark soliton. Quasi-one-dimensional traps are realistic prospects in the 
: relatively near future [H, and dark solitons are expected to emerge in them from generic violent collisions between 
' condensates [pMI- A controlled method for creating them by adiabatic state engineering with lasers has also recently 
, been proposed Jo] . And they are expected to play a crucial role in the eventual decay of superfluid currents in tight 
OsJ ■ toroidal traps which would be a valuable analog of the thin superconducting wires whose resistivity is one of 
the triumphs of non-equilibrium statistical mechanics Although dark solitons have been studied extensively in 
nonlinear optics optical fibres are spatially homogeneous on the relevant scale. In this Letter we extend or correct 
previous treatments of dark soliton motion in Bose condensates [p|,ppC|], by using multiple scale analysis to derive 
equations of motion for a dark soliton moving through a background condensate which changes slowly in both space 
' and time, and is subject to a generic slowly- varying potential (not necessarily harmonic). This powerful analytical 
^— I method may also be useful for other structures. 

The Gross-Pitaevski equation (GPE) governs the evolution of the c-number 'macroscopic wave function' tp{x,t) 
of a Bose-Einstein condensate. (This is of course a mean field approximation to the full quantum field theory; we 
• will consider dissipation from quasi-particle interactions very briefly below.) Incorporating a chemical potential by 
extracting a factor e"*''*/'', and then appropriately scaling the wave function, space, and time, one can write this 
equation in the convenient form 



> 
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^ ' We assume here a positive scattering length; and we do not restrict the normalization constant U = Jdx]^]"^, which 
is the number of particles rescaled by the strength of their mutual repulsion. Crucially, we assume a trap so thin 
that one can apply the GPE in one dimension. The approach to this limit from three dimensions has recently been 
discussed II 1^]. The essential requirement is that the transverse thickness of the trap be less than the healing length, 



to stabilize against buckling modes in the GPE. Making transverse confinement stronger than the temperature will 
' make even the quantum field theory effectively one dimensional. Experimental capability is already approaching both 
these limits. 

Eqn. ([^) in one dimension with constant V has been extensively studied in nonlinear optics and a solution 
with a localized structure has long been known: iIjus = tanh VI — V x. This time-independent solution is known as a 
dark soliton, because it describes a small dark spot in a light pulse; in our case this becomes a small 'bubble' of low 
condensate density in the dilute Bose gas. If tp{x) were restricted to be real, the dark soliton would be topologically 
stable, like other 'kink' solitons; but by taking -0 into the complex plane one can deform it into a configuration 
with constant density and phase, eliminating the 'bubble'. So unlike two-dimensional vortices, dark solitons are not 
topologically stable. Before considering their motion, therefore, one should first examine their stability; but in fact 
the two problems are closely connected, because the c omp lex deformations of ipos include the larger family of dark 
solitons moving with arbitrary (sub-critical) velocities |l3|] . These exact non-stationary solutions to (|l|), moving with 
constant velocity p= q, are 



V'GS = ip+ \Jvl-pP- tanh ^ - jp-{x - q[t)\ , (2) 



where = \f\--V is the Landau critical velocity (which the soliton cannot exceed). For p ^ we recover the 
motionless dark soliton at position x = q. Since for non-zero p the condensate density never vanishes, moving 
dark solitons are also called grey solitons. 
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For moving solitons the difference between maximum and minimum densities is v"^ — P^, and the phase shp across 
the sohton is tt + 2 arctan(p/-\/ v"^ — p^)- This means that in the hmit p —^ ±Vc, the sohton becomes identical with 
motionless condensate. Thus the soliton with maximum speed is the ground state; the energy of slower solitons is 
higher! In this sense one may say that dark solitons have negative kinetic energy. To be precise: The solutions given 
in ^ have fixed chemical potential (set to one), and their free energy G = E — U is 

G = Go + j{l-V~p')K (3) 

where E = ^ J dx [lip' + + 2V\^p\'^], and Go is the free energy of the ground state = Vl ^ V- (Alternatively we 
could change the chemical potential with p so as to keep the particle number constant; this slightly different family of 
grey soliton solutions with constant U has energy E = i?o + (4/3)(l — V—p'^)^^'^.) Thus, dark solitons are energetically 
as well as topologically unstable; but their instability is to acceleration, not to filling or collapse. Bogoliubov theory 
shows that acceleration is indeed their only instability, and as we discuss below, the 'anti-damping' time scale should 
be quite long. 

So apart from slow anti-damping, dark solitons in bulk behave as robust free particles, obeying q = p = 0. We now 
consider a dark soliton in a slowly varying medium, where we will be able to derive a more complicated equation of 
motion, if we interpret the 'slow variation' of V{x) as implying that there exists a length scale A which is both large 
compared to the soliton sca le and small compared to the trap scale. Precisely: there is a small dimensionless e, such 
that exp(— ^1 — V{q) ~ p'^K) << e for all phase space points {q,p) through which the soliton will actually pass, but 
V{x) = V{q) + V'{q){x — q) + C(e^) as long as |x — < A (with y((7)|A| being of order e). We will then examine 
an interval |a; — (/j < A around a grey soliton in a trap, the interval moving with the soliton, and smoothly patch 
this interval into a background condensate cloud in the hydrodynamic limit. Applying a simple form of multiple time 
scale analysis will then yield the equation of motion. This involved procedure ('boundary layer theory') is indeed 
necessary: merely treating V{x) as a perturbation is only valid if the potential is everywhere small, whereas we are 
interested in cases where, over large enough distances, it can change greatly. And ordinary perturbation theory will 
be valid only for a short time, but we are interested in large changes over longer times (such as the reflection of the 
soliton from a barrier). 

We begin with the simplest step of considering the background cloud. We will assume that the background cloud 
consists of condensate varying slowly on the healing length scale and its associated time scale, except possibly for 
small high frequency perturbations. For the dominant low-frequency component, we define ^p = ^/pe^^ for real p, 9, 
and stipulate that spatial and temporal derivatives of p and v = dx9 are of order e. We may therefore neglect y^"/ ^/p 
in the GPE to obtain the hydrodynamic equations 

dtp = -dxipdxO) = -dx{pv) 

dtO =1- p-V - . (4) 

We now patch our family of solitons into this background condensate: within |a; — gj < A we write ■0 = 
e'^+™(^-«)[?Ao + eV'i(a;- g,i)] for 

= i{p — v) + n tanh Hi{x — q) , (5) 

where 9{t) = [9{q — A,t) + 9{q + A, t)]/2, p and v are similarly defined, and n"^ = p — {p — v)"^ . Since p would be 
constant if e — *■ 0, we conclude that p is order e; in fact, p, v, and k may be taken as functions of the 'slow time' et. 

We then expand the Gross-Pitaevski equation to order e within |a; — <z| < A, keeping in mind that p, V' , k etc., are 
all 0(e). First using (0) to estabhsh 

dt9 



we find from the zeroth order terms q - 

[V\q)+ij](x- 

We will abbreviate Eqn. (^ as idt^i\x-q + 5(V'i,'!/'i) = J{x,et). (Note that it is a straightforward but very impor- 
tant step in obtaining (R) to distinguish 9t, which is, as usual, differentiation with respect to t with x fixed, from 



= \ll[^dt+pdx)9{x^t)\^^^^^ 

= pv + l- p-V{q)-v^/2 + 0{e^) , (6) 
p as always, plus the following at first order in e: 
. d 

- q)ipo + [p — V — i/i— K tanh hi(x — q)] 

Hi 

- i{p - vW, + - ~ (2|7^oP - P)^i] . (7) 
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differentiation with respect to t with x — q fixed: dtf{x — q,t)\x = dtf{x — q,t)\x-q — pdxf{x — q,t)\t for any function 

/■) p. . . 

We could then proceed to solve Eqn. (Uu using the Green's function for the homogeneous part. To construct this 

we would need all the independent solutions to the homogeneous equation; but in fact for our purpose we will require 
only the four independent solutions ui, U4 to the time-independent equation £{uj, u*) = 0. Distinguishing the fast 
and slow parts of -01 by defining ipi = (j){x — q, et) + xi^ ~ 9) we can use (Q) to show that the real parts of certain 
integrals are constrained to vanish: 

Re / dx [2iu*{x)dtx\x-q + d^lx'u* -u* x- 2i(p - v)u*x] 

J q-A ^ 

= Re^'^ dx (2u*{x)J{x, et) - dx[(j)'u* - u* 4> - 2i{p - v)u*(j)]j = . (8) 

Here the crucial final equality follows from the fact that the two sides of the preceding equation vary on different time 
scales, and so must separately equal zero. (Since the first line is linear in x, which must be fast, a non-zero constant 
is not allowed.) This is the great strength of the combination of boundary layer and multiple time scale analysis, that 
it allows us to obtain the motion of a short-scale defect in a long-scale background, by solving only time-independent 
equations. 

Eqn. (H) gives us four constraints, which since all four Uj{x) may be obtained explicitly, can be evaluated. In 
addition we require that our soliton il) match smoothly into the background flow as jx — g| — > A, and this introduces 
constraints from (^) as well. Together these constraints fix the hitherto unknown p, and also relate p,6,v aX x — q~ K 
to their values aX x = q + K. We illustrate the procedure with the simplest but most important constraint, the 
one involving ui{x) = sedc? k{x — q). Since ui(±A) and u[{±A) are exponentially negligible, we discard terms of 
this order in (^). We can then extend the limits of integration to infinity and shift the integration dummy variable 
X — g — !■ y, to obtain 

dy [K,{V'{q) + v)y tauh ny + p — v] sech^Ky = 2p + V'{q) — v — . (9) 

This is the equation of motion, accurate to C(e), for a dark soliton in an otherwise hydrodynamic condensate in an 
inhomogeneous potential. We will examine it in some simple limits, before discussing the conditions obtained from 
the other Uj, and from requiring (^) as |a; — qj — > A. 
With V = 0, Eqn. (^ imphcs 

q = -lv'iq). (10) 

In a harmonic trap, this implies oscillation of the soliton with frequency l/\/2 times that of the dipole mode of 
the condensate (the trap frequency) ||T^. This result can also be obtained for small oscillations by solving the 
Bogoliubov equations for a motionless soliton in a trap, using a simpler, time-independent version of the 'boundary 
layer' approach that led to (joj) |jl^. Wc have confirmed this frequency to rather more than the expected accuracy 
in numerical simulations fl^ ] of harmonic traps over a wide range of condensate densities and oscillation amplitudes; 
we have also confirmed that the center of mass is decoupled and oscillates at the trap frequency. Eqn. ( p^ ) also 
holds for arbitrary potentials, however, as long as they vary slowly on the healing length scale. We have therefore 
further confirmed the good accuracy of our equation of motion by solving Eqn. (|l|) numerically over a wide range 
of parameters and for various potentials; a generic example is shown in Fig. |l|. Since with lasers one can generate 
micro-wells or barriers in a trap, it should be possible to realize similar potentials experimentally. 

We now consider a stationary background flow, such as in an inhomogeneous toroidal trap holding a persistent 
current. In general the system is quite complicated; but in the limit where both the inhomogeneous potential V and 
the average kinetic energy Vq are small compared to the chemical potential, we have p = l — V, v = vo[l + V], which 
with dtv = implies the easily solvable equation 

q = V'{q)[voq-l]/2. (11) 

Despite the q term, Eqn. ( |Tl| ) is not dissipative: it may be derived variationally from the Lagrangian (2/tig)(l — 
wo9)[ln(l — voq) — 1] — V, and the energy q^ — L is conserved. 

A simple example of the generally still more complex case where p and v are time-dependent is a soliton moving in 
a harmonic trap of frequency fl in which the collective dipole mode has also been excited: 
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q = - — [q + Qcosn{t-to)], (12) 

where Q is the dipole amphtude. As required by the Ehrenfest theorem for a condensate in a harmonic trap, the rigid 
dipole oscillation of background and soliton together, q = Q cos il{t — to), is a solution to (p^. 

Since this Ehrenfest theorem states that the centre of mass of the condensate must oscillate at the trap frequency 
il, but ( |To| ) makes the small soliton 'bubble' oscillate at fl/y/2, it is clear that the background condensate must be 
perturbed by the soliton moving through it. This brings us back to the constraints we have not yet examined, which 
turn out to imply discontinuities of 0{e) in both p and v between x — q = ±A. These are in addition to the trivial 
discontinuities oc A due to background gradients. It is both convenient, and consistent with our ©(e), 'boundary layer' 
approach, to consider the entire interval |a; — q| < A to be pointlike as far as the background condensate is concerned; 
so, formally letting A ^ after obtaining all our results so far, the discontinuities across the soliton become abrupt. 
The requirement for them can then be expressed as delta function sources, ai x = q(t), which must be added to the 
hydrodynamic equations. The result can be shown to be 

dtp = —dx{pv) + 2kS{x — q) 

dtv = -dx{v^/2 + p + V)+ p-^S{x - q)[K{V' + v) + 2{p - v)k] . (13) 

In most cases indeed these delta function sources are unimportant, since the soliton couples only to the smooth part 
V = [v{q+,t) + v{q—,t)]/2, and the sources generate only discontinuities. The effect of these on v depends on the 



boundary conditions for the entire condensate, and solving (|9D and (|13| ) together to determine this effect is generally 
not much easier than numerically solving the GPE with the dark soliton. There are nevertheless some important 
points that can be learned from the source terms. For instance, they preserve the Ehrenfest theorem in a harmonic 
trap, as may be checked straightforwardly by evolving X = —2nq + Jdxxp under (^) and (p^). And because of its 
coupling to the background fluid, one can deduce that a dark soliton oscillating in a small well within a large sample 
of bulk condensate will generate sound waves, and so exhibit radiative anti-damping. Numerical integration of the 
GPE confirms this prediction: the soliton eventually escapes from the micro-well, the radiation ceasing as it enters 
the region of constant potential |l^ . 

In a finite trap, however, coupling to the background condensate modes does not provide dissipation. In this case 
dissipation can only come from corrections to mean field theory; in particular, from collisions with uncondensed 
atoms of the thermal cloud. A simple estimate of the anti-damping time scale is provided by the rate at which the 
soliton encounters particles, divided by the number 2k of particles 'in' the soliton (for the 'soliton mass'). At current 
experimental temperatures and densities, with 99% of the particles in the condensate, this time is on the order of one 
second; which agrees with the calculation in Ref. |jl^ of the dark soliton decay time. It is clear therefore that the 
instability of dark solitons is by no means fast enough to prevent their observation. 
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FIG. 1. Density for a dark soliton oscillating through a static Thomas-Fermi cloud at J da;|'!/'|^ = 300, with potential 
V = 0.1a;(a; — 2) + l.lsech'^a; shown in dots. Initially q = —2.18. Eqn. ( [lo| ) thus predicts the time between turning points 
(where — — Q aX the minimum) to be T/2 — 11.5; the error of about 4% is indeed 0{e). 
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